NSite = Nx/lc;
beta = 1000;
Source = 1 : Nx/4;
Target = Nx/4*2+1 : Nx/4*3;
% Source = 1 : Nx/8;
% Target = Nx/8*2+1 : Nx/8*3;

mu = D(NSite);
RhoExact = real(2d0 * EF * diag( 1./ (exp(beta*(D-mu))+1d0) ) * EF');
S = svd( RhoExact( Source, Target ) );
semilogy(S, 'b');
hold on

mu = D(2*NSite);
RhoExact = real(2d0 * EF * diag( 1./ (exp(beta*(D-mu))+1d0) ) * EF');
S = svd( RhoExact( Source, Target ) );
semilogy(S, 'r');


mu = D(4*NSite);
RhoExact = real(2d0 * EF * diag( 1./ (exp(beta*(D-mu))+1d0) ) * EF');
S = svd( RhoExact( Source, Target ) );
semilogy(S, 'k');

